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Abstract 

From the Schwarzschild metric we obtain the higher-order terms (up to 20-th or¬ 
der) for the deflection of light around a massive object using the Lindstedt-Poincare 
method to solve the equation of motion of a photon around the stellar object. Ad¬ 
ditionally, we obtain diagonal Pade approximants from the perturbation expansion, 
and we show how these are a better fit for the numerical data. Furthermore, we use 
these approximants in ray-tracing algorithms to model the bending of light around 
the massive object. 
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Introduction 


The General Theory of Relativity (GTR) is probably one of the most elegant theories ever 
performed. It was put forth by Albert Einstein in its current form in a 1916 publication, 
which expanded on his previous work of 1915 in Ei- This is summarized in 14 equations 
M- The Einstein field equations (ten equations written in tensor notation) 


G fiu — Rfiu 


2 ^9 i/ T Xq^i/ 


and the geodesic equations (4 equations) 


(i) 
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( 2 ) 


In Q G^ is the Einstein’s Tensor, which describes the curvature of space-time, R IJV is 
the Ricci tensor, and R is the Ricci scalar (the trace of the Ricci tensor), g^ is the metric 
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tensor that describes the deviation of the Pythagoras theorem in a curved space, is 
the stress-energy tensor describing the content of matter and energy, k = where 
c is the speed of light in vacuum and G is the gravitational constant. Finally, A is the 
cosmological constant introduced by Einstein in 1917 ME! that is a measure of the 
contribution to the energy density of the universe due to vacuum fluctuations. In equation 
([2J) s is the arc length satisfying the relation ds 2 = g^dx^dx" and Y^ a are the connection 
coefficients ( Christoffel symbols of the second kind), ad 1 is the position four-vector of the 
particle. We use Greek letters as n,u,a,etc for 0,1,2,3. We have adopted the Einstein 
summation convention in which we sum over repeated indices. Einstein’s equations ( [Tj) 
tells us that the curvature of a region of space-time is determined by the distribution of 
mass-energy of the same and they can be derived from the Einstein-Hilbert action [8, 9j: 



( 3 ) 


where 3ft represents a region of space-time, L F is the Lagrangian density due to the fields 
of matter and energy and g is the determinant of the metric tensor. 

One of the most relevant predictions of General Relativity is the gravitational de¬ 
flection of light. It was demonstrated during the solar eclipse of 1919 by two british 


expeditions [ID] • One of the expeditions was led by Arthur Eddington and was bound for 
the island of Principe in East Africa. The other one was led by Andrew Crommelin in the 


region of Sobral in Brazil. The light deflection can be measured taking a photograph of 
a star near the limb of the Sun, and then comparing it with another picture of the same 
star when the sun is not in the visual field. The observations are not easy. At present, 
Very Long Baseline Interferometry (VLBI) is used to measure the gravitational deflection 
of radio waves by the sun from observations of extragalactic radio sources HU. The result 
is very close to the value predicted by General Relativity m , which is Rec i 1.752 
seconds of arc (M@ and Rq represent the solar mass and radius, respectively). 

In the literature we can find calculations to second order of the deflection of light by a 
spherically symmetric body using Schwarzschild coordinates [12| [L3j, 14j [15] I 11 this paper 
using the Schwarzschild metric we obtain higher order corrections (up to 20-th order) 
for the gravitational deflection of light around a massive object like a star or a black 
hole using the Lindstedt-Poincare method to solve the equation of motion of a photon 
around the stellar body. Additionally, we obtain diagonal Pade approximants from the 
perturbation expansion, and we show how these are a better fit for the numerical data. 
We also use these approximants in ray-tracing algorithms to model the bending of light 
around the massive object. 



3 


1 Schwarzschild metric 

For a spherical symmetric space-time with a mass M in the center of the coordinate 
system, the invariant interval is mm 


(ds) 2 = 7 (cdt) 2 — 7 1 (dr) 2 — r 2 (dtt) 2 (4) 

where (dtt) 2 = ( d6 ) 2 + sin 2 9 (dtp) 2 , with coordinates x° = ct, x 1 = r, x 2 = 9 and x 3 = <p. 
7 = 1 — — where r s = is the Schwarzschild radius. 

The corresponding covariant metric tensor is given by 




7 0 0 0 

0 -7- 1 0 0 

0 0 -r 2 0 

0 0 0 —r 2 sin 2 9 


(5) 


Equation Q has two singularities. The first one is when r = r s (the Schwarzschild 
radius) which defines the horizon event of a black hole. This is a mathematical singularity 
that can be removed by a convenient coordinate transformation like the one introduced 
by Eddington in 1924 or Finkelstein inl958 m- 


t = t±— ln\— -II ( 6 ) 

c r a 

With this coordinate transformation the invariant interval reads: 

(ds) 2 = c 2 ^1-- j (dtj — ^1 H— (dr) 2 =f 2c didr — r 2 (dQ) 2 (7) 

The other singularity in r = 0 is a mathematical singularity. For a radius r < r s , 
all massless and massive test particles eventually reach the singularity at r = 0. Thus, 
neglecting quantum effects like Hawking radiation m Fl8j . any particle (even photons) 
that falls beyond this Schwarzschild radius will not escape the black hole. 
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2 Geodesic equation for a photon in a Schwarzschild 
metric 


The geodesic equation can be written in an alternative form using the Lagrangian 


L 



-9af} M 


dx a dx 13 

da da 


( 8 ) 


where a is a parameter of the trajectory of the particle, which is usually taken to be the 
proper time, r, or an affine parameter for massless particles like a photon. The resulting 
geodesic equation is: 


^ = \ (tV/afl) U a vP (9) 

where = d ^ L . 

Consider a photon traveling in the equatorial plane ( 6 = 7t/2) around a massive object. 
For a photon, dr = 0 and thus, we use an affine parameter, A, to describe the trajectory 
instead of the proper time, r. For the coordinates ct (fi = 0) and 4> (/r = 2) the geodesic 
equation ([9]) give us, respectively : 


d 2 dt 
dX 7C dX 


0 


and 


d 2 d0 
dX dX 


( 10 ) 


( 11 ) 


Both of these equations define the following constants along the trajectory of the 
photon around the massive object: 


^ E ' 


( 12 ) 


and 


2 d(j) 


r ~x = J 


(13) 


where E' has units of energy per unit mass and J of angular momentum per unit mass 
(when A has units of time). 
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The invariant interval for the Schwarzschild metric in the plane 9 = 7t/2 is. 


( ds) 2 = c 2 (drY = 7 c 2 (dt) 2 — 7 1 ( dr) 2 — r 2 (dp) 2 = 0 . 


(14) 


Using 7 x ~ the equation can be written in the form 


7 c 


2 (dt\ 2 ^ fdr\ 2 /70 \ 2 _ r2 /70 \ 2 = Q 


dX 


7 


d(fi I \d\ 


dX 


(15) 


Multiplying (15) by 7 , and inserting the definitions of E' and J we obtain: 


(E'f _ J 2 (dr\ 2 _ 7 J 2 
c 2 r 4 1 ) r 2 


= 0 


(16) 


This equation can be turned into an equation for U((p) = r( U, noting that 


dU 1 dr 

d(p r 2 d<fi 

so we arrive at the following equation for U((p): 


(17) 


7-7S1 


(18) 


By taking the derivative of equation (18) with respect to p, we get the following 
differential equation for U (0) 


s)(=S*»-«7" 


(19) 


The differential equation in (19) can be separated into two differential equations for 


17(0). The first one is the equation for a photon that travels directly into or out from the 
black hole: 


dU 

70 


= 0 


( 20 ) 


the other differential equation, applicable for trajectories in which 17(0) is not constant 
with respect to 0, is the following: 


d 2 U rr 3 tt2 

W +U= -2 r - U 


( 21 ) 


This equation can also be written in the following way, using the definition of the 
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Schwarzschild radius: 


d 2 U 


U = 


3GMU 2 


( 22 ) 


This is the equation for the trajectory of a massless particle that travels around a 
black hole in the equatorial plane. 


3 Differential equation for the trajectory of a photon 

In the previous section, we obtained a differential equation for a photon traveling around 


a rnasive object like a star or a black hole (see equation (22)). This equation has an exact 
constant solution, for the unstable circular orbit of a photon around the black hole: 


3 GM 


r r = 


(23) 


where r c is the radius of the so-called photon sphere [TB]. We note that the radius of 
the photon sphere can be expressed in terms of the Schwarzschild radius: 


3 r. 


r, = 


(24) 


The orbit described by a photon in the photon sphere is actually an unstable orbit , 
and a small perturbation in the orbit can lead either to the photon escaping the black 
hole or diving towards the event horizon pL 6 ] . 


Equation (22) is nonlinear, and is highly difficult to solve analytically. However, a 


perturbative solution of this equation can be readily obtained. Let’s first rewrite Equation 


( 22 ) in terms of r c : 


d 2 U 

d(p 2 


U = r r U 2 


(25) 


Consider the initial conditions shown in Figure [TJ The smallest value of the r—coordinate 
in the trajectory, r = b, is taken such that the photon escapes the black hole, b > r c . We 


will rewrite Equation (25) in terms of e = ^< 1 , which we will use as a non-dimensional 


small number for our following perturbative expansions. Note that by multiplying both 


sides of Equation (25) by b, and defining the non-dimensional trajectory parameter 


VM = 


r(<p) 


(26) 


Equation (25), with the inclusion of the term e = W then becomes a differential 


equation in V (0): 
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Figure 1: Trajectory of a photon outside the photon sphere. The initial conditions are 
taken such that r |^ =0 = b, and is called the impact parameter of the trajectory - the closest 
distance from the trajectory to the center of the black hole. We thus have, ^|</,=o = 0 
and jff |c/>=o = 0. The photon experiences a total angular deflection of 2a. 


where 0 < e 


d 2 V 

~d^ 


+ U 


eV 2 


y < 1, and with initial conditions given by 


1/(0 = 0 ) 



(0 = 0 ) = 0 


Under these conditions, U(0) is bounded such that 


(27) 


(28) 


\V(4>)\ < 1 


( 29 ) 





4 First-order solution for V{(j>) 


A first idea to obtain a solution of Equation (27) is to consider a V(<f>) as a power series 
in e: 


V{fr e) = E o ( 0 ) + eK( 0 ) + eV 2 (0) + ... 


Plugging the expansion (30) into Equation (27) results in the following: 


(30) 


'd 2 V 0 d 2 V l 2 d 2 E 2 


+ ... j + (Vq + ehi + e 2 V 2 + •••) — £ ^Vq + eV7 + e 2 !/^ + ...^ (31) 


We can group the powers of e in Equation (31): 


0 . d 2 Vo I t/ _ n 

e ■ - U 


^ + u = r„ ! 


f 2 : ^ + V 2 = 2V 0 V 1 


3 • + U = (V,) 2 + 2V 0 v 2 


(32) 

(33) 

(34) 

(35) 


Note that the initial conditions of V (</>), applied to the asymptotic expansion in Equa¬ 
tion (30), imply the following, by grouping powers of e: 


V„(0) = 1;S*(0)=0 (36) 

e*: n(0) = 0;!»(0) = 0;A:>l (37) 

From these differential equations and initial conditions, we can readily obtain Vo and 
V\ iteratively Q 


1 It is convenient to write the Vk((f>) in terms of polynomials in cos((f >) 
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Vq{ 4>) = cos{<f>) 


(38) 


V\((p) = ^ - ^cos(0) - ^COS 2 (0) 


Thus, we obtain an equation for V(4>), per Equation (30): 


V (0) = cos(0) + e 


J - 0os( <t>) - icos 2 (^) + 0(£ 2 ) 


(39) 


(40) 


According to the coordinate system shown in Figure 1, the photon goes through a 
total angular deflection of 2a. This corresponds to setting V(<f>) — 0 for both 0 = n/2 + a 
and 0 = — 7t/2 — a. From both of these conditions considering that a is very small, to 
first order in e we get: 


a = 


2e 

3 


(!-!) 


(41) 


The total deviation of the photon is then 


4e 4r c 4GM 

[} = 2a& — = —r = -r-=- 

3 3b be 2 


(42) 


For a light ray grazing the Sun’s limb b = Rq = 695510 km [T9] and we get the very 
well known value 


9 = 2 a 


4 GM e 
Rec 2 


= 1.7516 arcseconds 


(43) 


where M 0 = 1.9885 x 10 M kg is the Sun’ s Mass, and c = 2.99792458 x 10 b 
value of the speed of light in vacuum [19] . 


— is the 

S 


5 Towards a second-order solution for Q(e) 

We will now see how to obtain higher-order solutions for 9. The differential equation in 


(34) has the following solution: 


4 41 2 1 5 

^ 2 ( 0 ) = -g + ^ C0S (^) + g cos 2 (0) + —COS 3 (0) + —0sm(0) 


(44) 


However, the term in Equation (44) that goes as 0sin(0) grows without bound, and 


occurs because the right-handed side of Equation (34) contains terms proportional to the 


homogeneous solution of Equation (34): acos(4>) + b sin(0). When this happens, the 
solution contains terms that grow without bound, such as 0sin(0), called secular terms 
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PI Thus, if we naively include Equation ( |44| ) in V (0), our solution is no longer bounded. 
Thus, we have to eliminate any and all secular term that arises to arrive at a well-behaved 
solution for 1/(0). 

One method to do this, due to Lindstedt and Poincare, is by solving the differential 
equation in the following strained coordinate 


0 = 0 (l + cnie + cn 2 e 2 + ...) (45) 

Where the Uk are constants to be determined. In terms of this new strained coordinate 


0, Equation (27) becomes 


d 2 V 


(l +C o 1 e + cu 2 e 2 + ...) + V (4>) = eV 2 (4>) 

v ' d(p z 

We proceed in the previous way, and assume an asymptotic expansion on 1/(0): 


(46) 


V(4>- e) = V 0 (4>) + eU(0) + e 2 vS) + ... 


(47) 


Plugging the expansion(|47[l in Equation (461, we obtain: 


(1 


+ cuie + oj 2 e- + ... 


2 


d 2 V 0 d 2 Vx 2 d 2 V 2 

dcj) 2 d<t> 2 d<p 2 


+ 



(48) 


We can group the powers of e in Equation (48): 


e 


o . 


d 2 Vo 

d<fi 2 


+ Vo — 0 


(49) 


e 


i . 



V 0 2 - 2cni 


d 2 V 0 

d(f) 2 


(50) 


e 


2 . 


d 2 Vo 

dtp 2 


+ v 2 


2VqVi — ( 'ui\ 2 + 2oj 2 ) 


d 2 V 0 
d(p 2 


— 2 u 


d 2 V i 

^ dip 2 


(51) 


,3 . d 2 V 3 , t/ 

6 ' -W +V 3 


lb 2 + 2V 0 V 2 - (2 Ul u 2 + 2cn 3 )^ 


(cui 2 + 2cn 2 ) 


d 2 V ! 
dip 2 


— 2u 


d 2 V 2 
1 d<p 2 


(52) 
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With some care due to the definitions of the scaled variable and its derivative, we arrive 
at initial conditions for the 14(0) from the initial conditions of V(0): 


£ »: Vo (0 — 0 ) — 1 ; «( 0)=0 


(53) 


UW = 0) = 0; Ai(0) = 0;fc>l 


(54) 


Solving the differential equation (49) with initial conditions (53), we arrive at the 
zeroth-order contribution to 1/(0): 


V„(<j>) = cos(<j>) 


(55) 


Similarly, we can obtain 14(0) from Equation (50) subject to initial conditions (54): 


2 1 ~ 1 ~ „ 
W) = j - »«>s(<6) - -cos 2 ((f>) + ui<j>sin(<j>) 


(56) 


We note that a secular term has appeared for 14(0)- However, we use our freedom in 
the definition of uq to eliminate this secular term by setting 


U] — 0 

so that the final form of 14(0) is: 

14 ( 0 ) = ^ - ~cos(0) - ^cos 2 (0) 
Similarly we obtain for 14(0) : 


(57) 


(58) 


UOM = -g + jj cos (h + \co S \4>) + -^cos 3 (ij>) + 4j( 144w 2 + 60 )<j>sin(<j>) (59) 

To eliminate the secular term in 14(0), we set 

o; 2 = -^ (60) 

and obtain the well-behaved second-order term 

UM = -g + y°s(i) + J cos 2 (i ) + Tcos 3 (^) 


(61) 
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From all the solutions obtained so far, we can obtain the second-order correction to 0(e). 
Note that V (</>) is given by: 


+e 


(2 1 ~ 1 

V{4 >) = cos(4) + e (j - -cos(0) - -cos 2 (0) 

+ ^ cos (0) + + Y2 COs3 $)+) + °( e3 ) 


(62) 


We set up 0 = 7t/2 + a in Equation (62), such that V(tt/ 2 + a) = 0 and obtain: 


+e z 


■ /~x /2 1 . 1 . 2/ A 

— stn(a) + e I - + -sm(a) — -sm (a) I 
\ o o o / 

+ ^sin 2 (a) - -^sm 3 (d)+] + 0(e 3 ) = 0 
9 ob 9 12 / 


(63) 


We could truncate this Equation and solve the resultant cubic polynomial in sin(a). 
However, this method would not be easy to generalize, because we do not have a general 
formula for the roots of fifth-order polynomials and above, according to Galois theory 
EH- Also, an n —th order polynomial results in n different complex solutions, one of 
which we expect to have a leading term of order e, to obtain a better approximation of H, 
and we would need to check all the n different solutions for this. Additionally, we have 
to remember that so far this is an asymptotic expansion in e, and the truncation of the 
higher-order terms does not allow us to clearly see what the order of our estimate for ff(e) 
is. All of these problems are solved by assuming that sin(a) has the following expansion 
in e, with a leading term of order e 1 : 


sin(a) = exi + e 2 y 2 + e 3 y 3 + • ■ ■ (64) 

where the Xk are constants to be determined. Inserting this new expansion into Equa¬ 


tion (63) leads to the following algebraic Equation: 

'2 


(3- Xl ) e+ (-g + f -x,)t 2 + 0(J) = o 


(65) 


Then, we have to equal to zero the different powers of e in the last equation. Equating 
to zero the terms with e 1 we arrive at: 


2 

^ = 3 

and equating to zero the terms with e 2 we arrive at: 


X/2 = ~ 


2 

9 


( 66 ) 


(67) 
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Thus, sin(a ) is given by: 


sin(a) = - ^e 2 + 0(e 3 ) 


To obtain <3, we employ the Taylor series of arcsin(x ) around a; = 0: 


ar 


arcsin(x) = x + £4 + 0(a: 5 ) 
6 


( 68 ) 


(69) 


and obtain 


a =y-y 2+ °(' 3 ) 


(70) 


However, what we actually want is a. From the definition of the strained coordinate 


0 in (45), it is clear that: 


71 


2 T d \ 

2 ’ 01 1 + uqe + U7 2 ( 2 + 0(e 3 ) 

From the last equation, an using the Taylor expansion of — £^ 0 (—l) n a; n around 
x = 0, we obtain: 





e 2 + 0(e 3 ) 


From which we can obtain the total deflection angle, Vl = 2a 


(72) 





e 2 + 0(e 3 ) 


4 GM /15vr \ (GM\ 

6c 2 + ' 4 _ 4 J ' hr- ) 



(73) 


This result is in agreement with other work H2i [mm ees]. 


6 Higher order solutions for Q(e) 

The previous procedure can be automated to obtain higher-order expressions for Q. No¬ 
tably, all the solutions for the 14(0) are in the forms of (k + 1)—order polynomials of 
cos(0), and a secular term that is eliminated by choosing a suitable c Ok- The use of the 
expansion of sin(ci) in powers of e guarantees both the form of sin(a ) with a leading term 
of order e, and leads to algebraic equations for the Xk that are exceedingly easy to solve. 
Notably, getting a higher-order solution conserves the lower-order terms. Consider the 
formal Taylor expansion of O around e = 0: 


f! — fiqe + + ^3e 3 T ... 


(74) 
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A table of the coefficients K n of the series of Q in (74) can be found in Table [TJ These 
K n were found using the method of the previous sections, and obtaining the solutions up 
to V 20 (0). 



Exact value 

Numerical value 

Ki 

4 

3 

1.33333 


57t 4 

12 9 

0.864552 


122 5tt 

81 18 

0.633508 

K4 

385tt 130 

576 81 

0.494911 

K b 

7783 385ir 

2430 432 

0.403082 

Kq 

103565tt 21397 

62208 4374 

0.338319 

«7 

544045 85085tt 

61236 31104 

0.290571 

k 8 

6551545tt 133451 

1327104 8748 

0.254143 

Kg 

1094345069 1169918757T 

39680928 13436928 

0.225577 

^10 

22681108457T 1091492587 

143327232 22044960 

0.202655 

K\i 

33880841953 18553890355tt 

374134464 644972544 

0.183902 

^12 

3278312542505tt 627972527 

61917364224 3779136 

0.168300 

^13 

17954674772417 1514986498025tt 

58364976384 15479341056 

0.155132 

«14 

135335969751125?r 53937207017735 

743008370688 94281884928 

0.143875 

^15 

1532445398265737 1138317723327785?r 

1432594874880 3343537668096 

0.134145 

k 16 

1094325341294717675-rr 4027582104301883 

1711891286065152 2005632824832 

0.125654 

Kn 

2064610875963794827 128887453213429625tt 

545532128354304 106993205379072 

0.118179 

^18 

1263396148548501892925tt 2657173119021192719 

554652776685109248 371328591568896 

0.111548 

^19 

1085138496158025821251 399330245672667033725?r 

79959423384502272 92442129447518208 

0.105625 

^20 

218695963585074038928865tt 75186822805298075761 

26623333280885243904 2913501256925184 

0.100303 


Table 1: Coefficients K n of the series of in (74). For these coefficients, we report both 


the exact values and the numerical values with 6 significant figures. 


Clearly, as e —>■ 1, fl(e) —>■ 00 , because the photon starts going around the black hole 
as it starts closing in the photon sphere (b —y r c ). This means that 0(e) has a singularity 


at e = 1. The Taylor expansion of 0(e) around e = 0 that we found at Equation (74) does 
not return an estimate for the position of this singularity, because a polynomial does not 
have a singularity. However, we can obtain Pade approximants for O around e = 0, and 
these will return an estimate for the position of this singularity. 

As a small refresher on Pade approximants, we note their definition. A Pade approx- 
irnant of a function f(x) is a rational function f^ M \x) of the form: 


f [L\M\( s _ Q o + a i x + °2^ 2 + ••• + cilX L 

1 + b\x + b- 2 x 2 + ... + b M x M 


( 75 ) 
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where f{x) and f^ M \x) are equal in their first L + M + 1 derivatives around x = 0 
[22 > 23j- A diagonal Pade approximant f^(x) is a Pade approximant in which N = L = 
M. We can obtain the diagonal Pade approximants for up to N=10 with the Taylor series 
expansion for 0(e). For example, the OM(e) Pade approximant is given by: 


487T + (64 + 167T — 157T 2 ) e 

96 + (32 - 30vr)e ^ ^ 

The exact formulas for the Pade approximants of 0(e) are rather complicated because 
of the powers of 7r involved. Due to this, our work with Pade approximants will be 
purely numeric. All the Pade approximants 0^(e) have a singularity of order 1 at a 
position around e — 1. The position of this singularity, e s , is tabulated for the 10 Pade 
approximants in Table [2] 



N 

O 

1 

1.54222 

2 

1.21736 

3 

1.11036 

4 

1.06664 

5 

1.04532 

6 

1.03238 

7 

1.0245 

8 

1.01915 

9 

1.01537 

10 

1.01264 


Table 2: The position of the singularity near e = 1 for the Pade approximants, 0^1 (e). 


7 Numerical tests for Q(A) and its Pade approximants 


All the coefficients for the Taylor expansion of 0(e) were obtained around e = 0. We can 
test the correctness of the methods thus far used to obtain this function by comparing it 


to the results of numerical solutions of Equation (22). This is done with both truncated 


n —th order Taylor polynomials from O(e) and for the Pade approximants we obtain from 
this function, 0^(e). This comparisons are shown in Figures [ 2 ] and [ 3 ] 

We see from Figures [2] and [3] that the Pade approximants are much faster at converging 
into the actual form of 0(e). The convergence of the Pade approximants is such that for 
e = 0.99, the IV = 10 diagonal Pade approximant is within 3% of the corresponding 
numerical value. This is mainly due to the fact that Pade approximants are better in 
approximating functions that have singularities [23] . Once we know that the O(e) behave 
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0.0 0.2 0.4 0.6 0.8 1.0 

e 

Figure 2: Numerical points obtained for f2(e) compared to the truncated n —th order Tay¬ 
lor polynomials of O(e), up to 20—th order. With increasing value of n, the polynomials 
take larger values. 


correctly, we can use the O(e) to simulate the bending of light around a black hole. A 
simple first-order ray tracing algorithm that does this for the different approximations of 
O(e) we have found is shown in Appendix. 


8 Conclusions 


One of the most important predictions of the General Theory of Relativity is undoubtedly 
the bending of light around a massive object like a star, a black hole or even a galaxy 
(in this case it can generate a gravitational lens) . In this paper using the Schwarzschild 
metric we have obtained higher order corrections for the gravitational deflection of light 
around said objects using the Lindstedt-Poincare method to solve the equation of motion 
of a photon around the stellar body. We have successfully obtained an expression for (e) 
, the angular deflection experienced by a photon traveling around the massive object. We 
have assumed that the parameter e was small, and we were able to obtain the coefficients 
K n of the series of Q (e) up to V 2 o( 0 ) (the non- dimensional trajectory parameter, see 
equation (26)). The results are given in Table [T} Additionally, we have obtained diagonal 
Pade approximants from the perturbation expansion, and we have shown how these are 
a better fit for the numerical data. The best approximation for (e) we obtained was 
consistent with the numerical data even for an e ~ 0.99. In this case, the IV = 10 diagonal 
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Figure 3: Numerical points obtained for il(e) compared to the truncated N —th diagonal 
Pade approximants of h2(e), up to N = 10. With increasing value of N, the Pade approx- 
imants take larger values. For e = 0.99, the N — 10 Pade approximant is within 3% of 
the numerical value of 0(e). 


Pade approximant is within 3% of the corresponding numerical value. We were able to 
use this estimate for fl (e) in ray-tracing algorithms to model the bending of light around 
the massive object. 








Appendix 

A Ray tracing using Q(e) 
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Consider an observer A immersed in a background distribution of far away light sources. 
This observer can obtain the angular position of every object in the sky, and determine the 
intensity of light that comes from every point in the sky, Ia(9, 0), in spherical coordinates. 
Now, imagine another point in space, B , far enough from the observer A such that the 
intensity of light that comes from every point in the sky, according to an observer in point 
B , is also given by the distribution found by observer A: Ib(9, 0) = Ia(9, </>)• If we place 
a black hole at point B, then light coming from the faraway sources will bend around the 
black hole such that the original observer will see a different distribution of light around 
the black hole. In this condition, the observer will be able to note that the black hole 
effectively subtends a solid angle in the sky - region in the sky devoid of any light due to 
the black hole. One half of the angle subtended by the black hole will effectively give the 
’’angular radius” of the black hole, as seen by the observer, r B H- 

If we consider that the black region of the sky due to the black hole is due to the radius 
of the photon sphere, r c , instead of the the Schwarzschild radius, r s [^J and if we choose 
the coordinate system such that the black hole is at the positive x-axis, 9 = 7 t /2 and 
0 = 0, then, by definition of tbh , the new distribution of light measured by the observer 
will obey (for small enough vbh)- 

1(6 , 0) = 0 ; (9 - vr/ 2) 2 + 0 2 < (r BH ) 2 (A.l) 

For other values of ( 9 , 0), the observer sees light distribution shifted by the f2(e), where 
e is given by: 


Tc _ r BH 

b ~ ^(0 _ tt/2) 2 + 0 2 


(A. 2 ) 


for 9 a 7 r/2. We can use a further simplification of this latter equation, and use the 
coordinates (9 x ,9 y ) defined by 9 X = (ft, 9 y = 9 — 7 r/ 2 . For small values of 9 X and 9 y , say, 
in the order of milliradians, we can write: 


I(9 x ,9 y )=0;9 2 x + 9 2 y <(r BH ) 2 (A.3) 

2 One can convince himself of this by considering the light from faraway objects that grazes the black 
hole at a distance given by r = b. This trajectory of this light is bended by the black hole for b > r c . 
However, if b < r c , the light will not escape and effectively no light coming from faraway objects will 
seem to originate from r < r c , which becomes an effective radius for the black hole, according to observer 
A in these conditions. In the case that mass enters the black hole, and emits light from an r that obeys 
r s < r < r c , the light can escape the black hole, and is severely red-shifted. However, we are here 
considering a black hole with no light sources between r s < r < r c . 
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and 


^ _ r BH 

where the analogue with Cartesian coordinates is evident. This coordinate system is 
shown in Figure [I] for a black hole that subtends 47T x ICC 6 steradians, such that r BB = 2 
mrad. 

Gy / mRad 
10 



-5 


-10 

Figure 4: A black hole with r BB = 2 mrad in the center of the (6 X , 6 y ) coordinate system. 

This coordinate choice allows one to define the distribution of intensities that observer 
A sees to be (disregarding some attenuation factors): 



i{e x ,Q y ) = iA\e x -si{6) 




,6 y - fl(e) 




]0 2 x + 0 2 y < (r BH f 


(A.5) 


where we have used /a(^x, 9 y ), the angular distribution of intensities seen by observer 
A without the black hole present, and using the coordinates (9 X , 6 y ). We can see the effect 
of applying equation (A.5) by using the lA(9 xl 0 y ) defined from Figure [5] 

To model the deflection of light with the distribution in Figure[5j we use Equation (A.5) 
with f2(e) approximated as a truncated first-degree Taylor polynomial, and as diagonal 
Pade approximants with N = 2 and N = 10. The resulting images can be found in 
Figures [6] to [8] We use a black hole with r BH = 10 mrads. 
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Figure 5: 600 x 600 image corresponding to the intensity due to background light sources, 
Ia( 9 x , 9 y ), without a black hole present. Each pixel corresponds to 1 mrad. The big star, 
in white, has a radius of 50 mrad. The small stars, in gray, have a radius of 3 mrad. The 
star is at the center of the coordinate system, (9 x ,9 y ) = (0,0). 

The most notable difference between Figures [6] and [7] is the position of the white ring 
around the black hole, corresponding to the gravitational lensing of the big, white star at 
the black hole position 9 y = 0. When using a better approximation of fi(e), this ring has 
greater inner and outer radii, and is thinner. In Figure [8j there are 8 white pixels around 
9\ + 6y = 10 millirads, corresponding to a second ring of light due to the big, white star. 
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(C) (°) 

Figure 6: Background of Figure [ 5 ] warped by a black hole at (A) 9 y = 300 rnrad, ( B ) 
6 y = 200 mrad, (C) 9 y = 100 mrad, and (D) 9 y — 0 mrad. We make use of the truncated 
first-order Taylor polynomial of 0(e). 
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(C) (D) 


Figure 7: Background of Figure [ 5 ] warped by a black hole at (A) 9 y = 300 rnrad, (B) 
6 y = 200 mrad, ( C ) 9 y = 100 mrad, and (D) 9 y = 0 mrad. We make use of the diagonal 
N = 2 Pade approximant of 0(e). 
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(C) (D) 


Figure 8: Background of Figure [ 5 ] warped by a black hole at (A) 9 y = 300 mrad, ( B ) 
9 y = 200 mrad, (C) 9 y = 100 mrad, and (D) 9 y — 0 mrad. We make use of the diagonal 
N = 10 Pade approximant of 0(e). 
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